In [14]:
%pylab notebook
from charistools.convertors import Modice2Hypsometry
from charistools.readers import read_tile
from charistools.modelEnv import ModelEnv
import glob
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import re
pd.options.display.max_rows = 200


Populating the interactive namespace from numpy and matplotlib

In [8]:
%cd /Users/brodzik/projects/CHARIS/basins/major_basins_from_GRDC/MODIStiles
%ls


/Users/brodzik/projects/CHARIS/basins/major_basins_from_GRDC/MODIStiles
00notes.txt
00notes.txt~
AM_AmuDarya_at_Chatly_h22v04.tif
AM_AmuDarya_at_Chatly_h22v05.tif
AM_AmuDarya_at_Chatly_h23v04.tif
AM_AmuDarya_at_Chatly_h23v05.tif
AM_AmuDarya_at_Chatly_noholes_h22v04.tif
AM_AmuDarya_at_Chatly_noholes_h22v05.tif
AM_AmuDarya_at_Chatly_noholes_h23v04.tif
AM_AmuDarya_at_Chatly_noholes_h23v05.tif
BR_Bramaputra_at_Bahadurabad_h25v05.tif
BR_Bramaputra_at_Bahadurabad_h25v06.tif
BR_Bramaputra_at_Bahadurabad_h26v05.tif
BR_Bramaputra_at_Bahadurabad_h26v06.tif
BR_Bramaputra_at_Bahadurabad_noholes_h25v05.tif
BR_Bramaputra_at_Bahadurabad_noholes_h25v06.tif
BR_Bramaputra_at_Bahadurabad_noholes_h26v05.tif
BR_Bramaputra_at_Bahadurabad_noholes_h26v06.tif
GA_Ganges_at_Paksey_h24v05.tif
GA_Ganges_at_Paksey_h24v06.tif
GA_Ganges_at_Paksey_h25v05.tif
GA_Ganges_at_Paksey_h25v06.tif
GA_Ganges_at_Paksey_h26v06.tif
IN_Indus_at_Kotri_h23v05.tif
IN_Indus_at_Kotri_h23v06.tif
IN_Indus_at_Kotri_h24v05.tif
IN_Indus_at_Kotri_h24v06.tif
IN_Indus_at_Kotri_h25v05.tif
IN_Indus_at_Kotri_nolobe_gcs_noholes_h23v05.tif
IN_Indus_at_Kotri_nolobe_gcs_noholes_h23v06.tif
IN_Indus_at_Kotri_nolobe_gcs_noholes_h24v05.tif
IN_Indus_at_Kotri_nolobe_gcs_noholes_h24v06.tif
IN_Indus_at_Kotri_nolobe_gcs_noholes_h25v05.tif
IN_Indus_at_Kotri_nolobe_h23v05.tif
IN_Indus_at_Kotri_nolobe_h23v06.tif
IN_Indus_at_Kotri_nolobe_h24v05.tif
IN_Indus_at_Kotri_nolobe_h24v06.tif
IN_Indus_at_Kotri_nolobe_h25v05.tif
IN_Indus_at_Kotri_revised_h23v05.tif
IN_Indus_at_Kotri_revised_h23v06.tif
IN_Indus_at_Kotri_revised_h24v05.tif
IN_Indus_at_Kotri_revised_h24v06.tif
IN_Indus_at_Kotri_revised_h25v05.tif
SY_SyrDarya_at_TyumenAryk_h22v04.tif
SY_SyrDarya_at_TyumenAryk_h23v04.tif
SY_SyrDarya_at_TyumenAryk_h23v05.tif
basin_codes_h22v04.tif
basin_codes_h22v05.tif
basin_codes_h23v04.tif
basin_codes_h23v05.tif
basin_codes_h23v06.tif
basin_codes_h24v05.tif
basin_codes_h24v06.tif
basin_codes_h25v05.tif
basin_codes_h25v06.tif
basin_codes_h26v05.tif
basin_codes_h26v06.tif
empty_tiles/
indus_nolobe_qgis_project.qgs
indus_nolobe_qgis_project.qgs~
indus_revised_qgis_project.qgs
indus_revised_qgis_project.qgs~
mk_masks_IN_Indus_at_Kotri_new.sh*
tmp/

In [3]:
majorBasinIDs = ["SY_SyrDarya_at_TyumenAryk"]
#                 "AM_AmuDarya_at_Chatly",
#                 "AM_AmuDarya_at_Chatly_noholes",
#                 "IN_Indus_at_Kotri",
#                 "IN_Indus_at_Kotri_nolobe",
#                 "IN_Indus_at_Kotri_nolobe_gcs_noholes",
#                 "GA_Ganges_at_Paksey",
#                 "BR_Bramaputra_at_Bahadurabad",
#                 "BR_Bramaputra_at_Bahadurabad_noholes"]

In [ ]:


In [ ]:
def display_basin_mask(file):
    data = read_tile(filename=file)
    fig, ax = plt.subplots()
    print(np.amin(data), np.amax(data))
    ax.imshow(data, cmap='Greys', vmin=np.amin(data), vmax=np.amax(data), interpolation='None')
    ax.set_title(file)
    plt.axis('off')

In [5]:
for i in np.arange(len(majorBasinIDs)):
    
    list = glob.glob("%s_h[0-9][0-9]v*.tif" % majorBasinIDs[i])

    total_area = 0.
    for file in sort(list):
        print("file=", file)
        data = read_tile(filename=file)
        w2 = data[data == 1]
        print("shape=", w2.shape)
        area = w2.shape[0] * 463.312717 * 463.312717 / (1000. * 1000.)
        print("area= %f" % area)
        total_area = total_area + area
    
    print("%s total area= %f" % (majorBasinIDs[i], total_area))


('file=', 'SY_SyrDarya_at_TyumenAryk_h22v04.tif')
('shape=', (79258,))
area= 17013.417163
('file=', 'SY_SyrDarya_at_TyumenAryk_h23v04.tif')
('shape=', (977815,))
area= 209896.471057
('file=', 'SY_SyrDarya_at_TyumenAryk_h23v05.tif')
('shape=', (103223,))
area= 22157.712279
SY_SyrDarya_at_TyumenAryk total area= 249067.600499

Sanity check on MODICE area in major basin


In [15]:
configFile = '/Users/brodzik/ipython_notebooks/charis/calibration_modelEnv_config.ini'
myEnv = ModelEnv(tileConfigFile=configFile)
myEnv.tileConfig['input']['fixed']['basin_mask']


Out[15]:
{'dir': '%MODEL_TOP_DIR%/basins/major_basins_from_GRDC/MODIStiles/', 'pattern': '%DRAINAGEID%_%TILEID%.tif', 'id': 'CHARIS'}

In [16]:
from 
hyps = Modice2Hypsometry(majorBasinIDs[0],
                         myEnv,
                        contour_m=100,
                        modice_nstrikes=2,
                      verbose=True)


> charistools.modelEnv: fixed_filename is /Users/brodzik/projects/CHARIS/elevation_data/SRTMGL3_version2_SIN/CHARIS_DEM.v2.0.h22v04.tif; file_exists=True
> charistools.modelEnv: fixed_filename is /Users/brodzik/projects/CHARIS/basins/major_basins_from_GRDC/MODIStiles/SY_SyrDarya_at_TyumenAryk_h22v04.tif; file_exists=True
> charistools.modelEnv: fixed_filename is /Users/brodzik/projects/CHARIS/snow_cover/modice.v0.4/min05yr_nc/MODICE.v0.4.h22v04.2strike.min05yr.mask.nc; file_exists=True
> charistools.modelEnv: fixed_filename is /Users/brodzik/projects/CHARIS/elevation_data/SRTMGL3_version2_SIN/CHARIS_DEM.v2.0.h23v04.tif; file_exists=True
> charistools.modelEnv: fixed_filename is /Users/brodzik/projects/CHARIS/basins/major_basins_from_GRDC/MODIStiles/SY_SyrDarya_at_TyumenAryk_h23v04.tif; file_exists=True
> charistools.modelEnv: fixed_filename is /Users/brodzik/projects/CHARIS/snow_cover/modice.v0.4/min05yr_nc/MODICE.v0.4.h23v04.2strike.min05yr.mask.nc; file_exists=True
> charistools.modelEnv: fixed_filename is /Users/brodzik/projects/CHARIS/elevation_data/SRTMGL3_version2_SIN/CHARIS_DEM.v2.0.h23v05.tif; file_exists=True
> charistools.modelEnv: fixed_filename is /Users/brodzik/projects/CHARIS/basins/major_basins_from_GRDC/MODIStiles/SY_SyrDarya_at_TyumenAryk_h23v05.tif; file_exists=True
> charistools.modelEnv: fixed_filename is /Users/brodzik/projects/CHARIS/snow_cover/modice.v0.4/min05yr_nc/MODICE.v0.4.h23v05.2strike.min05yr.mask.nc; file_exists=True

In [17]:
hyps.data


Out[17]:
100.0 200.0 300.0 400.0 500.0 600.0 700.0 800.0 900.0 1000.0 ... 4600.0 4700.0 4800.0 4900.0 5000.0 5100.0 5200.0 5300.0 5400.0 5500.0
Date
NoDate 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 102.392187 65.041578 31.984142 14.382131 8.371688 2.790563 1.717269 0.858635 0.429317 0.214659

1 rows × 55 columns


In [ ]:
list = glob.glob("%s_h[0-9][0-9]v*.tif" % majorBasinIDs[1])

for file in sort(list):
    display_basin_mask(file)

#plt.close('all')

In [ ]:
df = pd.DataFrame(list, columns = ["File"])
df["basename"] = [re.search('(.+_OBJECTID[0-9]{1,3})', file).group(1) for file in df["File"]]
df["tileID"] = [re.search('(h[0-9]{2}v[0-9]{2})', file).group(1) for file in df["File"]]
df

In [ ]:
print(set(df["tileID"]))

In [ ]:
# Count the number of tiles per basin
basins = df.groupby(['basename'], as_index=False).count().sort_values(by=["File", "basename"])
basins.drop('tileID', axis=1, inplace=True)
basins.columns = ['basename', 'num_tiles']
basins

In [ ]:
# Get the basin prefix (basically the 2-letter ID) for output files
basin_prefix = re.search('(..)_', basins["basename"][0]).group(1)
for num_tiles in set(basins["num_tiles"]):
    out_filename = "%s_%dtile_basins.txt" % (basin_prefix, num_tiles)
    print("out file: ", out_filename)
    out_fh = open(out_filename, 'w')
    for subbasin in basins[basins["num_tiles"] == num_tiles]["basename"].values:
        # Get the list of tileIDs for this basin
        out_fh.write("%s\n" % subbasin)
    out_fh.close()

In [ ]:
%cat GA_4tile_basins.txt

In [ ]:
%cat SY_3tile_basins.txt

In [ ]:
%cat SY_2tile_basins.txt

In [ ]:
%cat SY_1tile_basins.txt

In [ ]:
# Special processing if we need to separate out basins that intersect only certain tiles
# Get the basin prefix (basically the 2-letter ID) for output files
basin_prefix = re.search('(..)_', basins["basename"][0]).group(1)
for num_tiles in set(basins["num_tiles"]):
    wh23_4_filename = "%s_%dtile_basins_with_h23-4.txt" % (basin_prefix, num_tiles)
    other_filename = "%s_%dtile_basins_other.txt" % (basin_prefix, num_tiles)
    print("next file: ", wh23_4_filename)
    print("other file: ", other_filename)
    h23_fh = open(wh23_4_filename, 'w')
    other_fh = open(other_filename, 'w')
    for subbasin in basins[basins["num_tiles"] == num_tiles]["basename"].values:
        # Get the list of tileIDs for this basin
        # Skip this basin if tileIDs include anything other than: [h23v05 h24v05]
        h23_4_basin = True
        for tileID in df.loc[df['basename'] == subbasin]["tileID"].values:
            if re.search('(h2[34]v05)', tileID):
                print("have this one")
            else:
                print("Don't have tileID=%s, need to skip basin=%s" % (tileID, subbasin))
                h23_4_basin=False
        if h23_4_basin:
            print("H23-4 BASIN: %s: %s" % (subbasin, ' '.join(df.loc[df['basename'] == subbasin]["tileID"].values)))
            h23_fh.write("%s\n" % subbasin)
        else:
            print("OTHER BASIN: %s: %s" % (subbasin, ' '.join(df.loc[df['basename'] == subbasin]["tileID"].values)))
            other_fh.write("%s\n" % subbasin)
    h23_fh.close()
    other_fh.close()

In [ ]:
%cat IN_2tile_basins_with_h23-4.txt

In [ ]:
%cat IN_2tile_basins_other.txt

In [ ]:
%cat IN_1tile_basins_with_h23-4.txt

In [ ]:
%cat IN_1tile_basins_other.txt

In [ ]: